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Cooperation and self-organized criticality are two main keywords in current studies of evolution. 
We propose a generalized Bak-Sneppen model and provide a natural mechanism which accounts for 
both phenomena simultaneously. We use the prisoner's dilemma games to mimic the interactions 
among the members in the population. Each member is identified by its cooperation probability, 
and its fitness is given by the payoffs from neighbors. The least fit member with the minimum 
payoff is replaced by a new member with a random cooperation probability. When the neighbors of 
the least fit one are also replaced with a non-zero probability, a strong cooperation emerges. The 
Bak-Sneppen process builds a self-organized structure so that the cooperation can emerge even in 
the parameter region where a uniform or random population decreases the number of cooperators. 
The emergence of cooperation is due to the same dynamical correlation that leads to self-organized 
criticality in replacement activities. 

PACS numbers: PACS Numbers: 02.50.Le, 87.23.-n, 87.23. Kg 



I. INTRODUCTION 

A fundamental question in the theory of evolution has 
been how cooperation can emerge between selfish mem- 
bers PUT] . Another question is why evolution takes place 
in terms of intermittent bursts of activities, which are 
the characteristics of dynamical systems in a 'critical' 
state 019]. Here, we propose a generalized Bak-Sneppen 
(BS) model [TU], which may solve the above two puz- 
zles simultaneously. We take an approach of evolutionary 
game theory and use the prisoner's dilemma (PD) games 
to mimic the interactions among members. Each mem- 
ber is identified by its stochastic strategy, specified by 
its (history independent) cooperation probability (CP). 
Here, a 'member' can represent an individual in a species, 
an agent in an economical system or a species in an eco- 
logical system. The fitness of a member is given by the 
payoffs of the games with its neighbors. We then ap- 
ply BS dynamics and replace the least fit member and 
its neighbors by new members with random CPs. The 
neighbors of the non-cooperator are likely to vanish due 
to its low payoff, but the non-cooperator itself can also 
be removed through the BS mechanism. As the non- 
cooperators disappear, the overall CP increases, and a 
new comer (with a random CP) will have a lower CP than 
the increased average. Therefore, the new comer tends to 
cause its neighbor to be the least fit, and the replacement 
activity likely occurs at or near the new comer's site. This 
invokes the spatio-temporal correlation between the least 
fit sites and can explain why replacements are episodic 
as well as how cooperation emerges. 

Evolutionary game theory has been one of the most 
powerful tools in studying the dynamics of evolution pQ . 
However, a simple straightforward application of game 
theory cannot explain the strong cooperation between 
"selfish" replicators observed in nature and society. For 
the evolution to construct a new, upper level of organi- 
zation, cooperation amongst the majority of the popula- 
tion is needed. However, the game theoretical descrip- 



tion of interactions between members usually leads to 
defections as evolutionarily stable strategies. Natural se- 
lection, which has been a fundamental principle of evo- 
lution, prefers the species that beat off the others and 
oppose cooperation. 



There have been numerous studies looking for natural 
mechanisms for the evolution of cooperation among com- 
petitive members HTHT^] . Recently, Nowak presented 
a state-of-art review on the evolution of cooperation and 
discussed five known mechanisms: kin selection, direct 
reciprocity indirect reciprocity, network reciprocity, and 
group selection [14] . Extensive studies provide the ex- 
act conditions for the emergence of cooperation for each 
of the five mechanisms. However, such conditions do not 
seem to be general enough to explain the cooperative phe- 
nomena observed everywhere. For example, for network 
reciprocity the benefit-to-cost ratio of a cooperative be- 
havior should be larger than the average degree 0, but 
this seems to be a rather strong assumption because the 
degrees are quite large in most cases in real population 
structures. Also, there have been a great deal of studies 
on self-organized criticality in game theory |15H19j . but 
their dynamics leading to the critical states are not di- 
rectly connected to the emergence of cooperation. Here, 
we consider an evolutionary game on networks and show 
that cooperation can emerge when the benefit to cost ra- 
tio is larger than just 1 if we use the BS process. When 
cooperators interact with defectors, they tend to disap- 
pear, giving rise to an assortment of cooperators |20) . 
Furthermore, this behavior emerges in the long run even 
with a small "chain-death" rate, w, where the number 
of neighbors that get replaced is less than one. For a 
uniform or random arrangement of cooperators and de- 
fectors, more cooperators than defectors disappear for 
small ui, but in the long run, the BS process builds a self- 
organized structure so that the number of cooperators in 
the population increases. 
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II. MODEL 

An influential model aimed to mimic the interactions 
between competitive members in a population is the PD 
game. It is one of the matrix games between two play- 
ers who have two possible decisions, cooperation (C) or 
defection (D). We consider a case in which the payoffs 
are calculated by the cost c and the benefit b of a co- 
operative behavior. If one player defects while the other 
cooperates, the defector receives benefit b without any 
cost whereas the cooperator pay cost c and its payoff be- 
comes — c. For mutual cooperation, both get benefit 6, 
but pay cost c, and their payoffs become b — c while the 
payoffs for mutual defection are 0. When we add c to 
all elements so that payoff can be directly interpreted as 
(non-negative) fitness, the payoff matrix becomes 

C D 
C ( b 0\ 
D \b+l 1 J ' 

where we set c = I without loss of generality. With con- 
ventional competition processes, the matrix game shown 
above does not, in general, predict the evolution of coop- 
eration. The birth-death process always predicts an evo- 
lution of defection. Cooperation can emerge for death- 
birth or imitation processes in a structured population, 
but only with a (unrealistic) large value of the benefit- 
to-cost ratio b for real populations [2J. 

Here, we consider the PD game interaction, but intro- 
duce the BS mechanism [10. as the competition process, 
and assume that the least fit member and its neighbors 
are prone to disappear. Each member is characterized 
by its strategy that determines when to choose the 'de- 
cisions' C or D. We consider the history- independent 
stochastic strategies, and the phcnotypc of a member, 
say the ith member, is represented by its CP Cj. The 
history independent pure (deterministic) strategies, the 
"always C" and the "always D" , correspond to the limits 
of Ci = 1 and a = 0, respectively. The fitness of a mem- 
ber is given by the sum of payoffs from its neighbors, and 
the member dies out if its total payoff is the minimum. 
The died-out site is occupied by a new member with a 
new CP, which is drawn randomly from to 1. Neighbors 
of the least fit site may also be harmed in the process of 
establishing the steady interaction with the new comer. 
Hence, we replace the neighbors of the least fit site by 
new members with the "chain-death" probability uj > 0. 



III. METHODS AND RESULTS 

We study the strategy evolution of a simple structured 
population from the initial state of random strategies. 
Initially, members in the population have cooperation 
probabilities that are drawn randomly from the uniform 
distribution of the interval [0 1]. They play PD games 
with their nearest neighbors. We assume that each mem- 
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FIG. 1: Real space configurations of (a) the CP Cj and (b) 
the RF fi for t € [0, 6000] with ui = 1 and N = 64. They are 
represented by colors, red for and green for 1, as indicated 
by the top panel. The black dots in (b) represent the least fit 
sites. 

ber plays sufficiently many games prior to the reproduc- 
tion process and use the payoff expectation value as its 
fitness. The least fit member with the minimum payoff 
expectation is replaced by a new member with a new 
random CP. In addition to the least fit member, the 
neighbors of the least fit member are also replaced by 
new members with the probability uj. Then, we recal- 
culate the payoff expectations, and replacements occur 
at the new least fit member and its neighbors. We con- 
tinue these processes until the system reaches a steady 
state and calculate the statistical properties of the popu- 
lation, such as the mean cooperation probability, fitness 
distribution, avalanche size (defined later) distribution 
and etc. 

For simplicity, we present our model and results in 
a one-dimensional (ID) structure, but our main results 
hold in other population structures. Initially (t = 0), we 
assign a random CP, c,(0), to the site i for i = 1, . . . , N. 
Then, we calculate the payoff expectation, /i(0), at time 
t = 0, 

fi(0) = b [ci_i(0) + Ci+i(0)] + 2 [1 - a(0)] , (1) 

of site i with a periodic boundary condition and find the 
minimum payoff site, mo- Except this minimum site, mo 
and its neighbors, m ± 1, the CPs are not changed at 
t = 1, so we set Cj(l) = c;(0) unless i — ttiq or mo ± 
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FIG. 2: Time dependence of MCP, G, for (a,c) ui = 0.2 and 
(b,d) u = 1.0 for systems with N = 32, 64, 128, and 256. In 
(a) and (b), the overall behaviors of MCP are shown while 
the initial transient characteristics are shown in (c) and (d). 
For uj = 0.2, MCP decreases first and then increases while it 
monotonically increases for u = 1. 



1. The CP at the m site, c m „(l), is given by a new 
random number between and 1. For its neighbor sites, 
Cm ±i(l) is given by a new independent random number 
with the probability uj, but remains as c mo ±i(0) with the 
probability 1 — uj. Now, we recalculate the payoff of 
Eq. ([!]) with Cfe(l) instead of Cfc(0). We find the new 
minimum payoff site, mi, of t = 1 and apply the same 
replacement dynamics to get t = 2 configurations and so 
on. 

Figure 1 shows typical real space configurations of CP, 
d, and the reduced fitness (RF), /, = / i /(26 + 2) € [0 1]. 
We show the configurations for initial 6000 time steps of 
a TV = 64 system with b = 1.5 and u> = 1. Both the 
CP and the RF are represented by colors, by red and 
1 by green. The least fit sites (black sites in (b)) and 
their two neighbors are where the replacement activity 
occurs. Comparison between the configurations in (a) 
and their equivalents in (b) reveals that the least fit sites 
are located where their neighbors are less cooperative 
[relatively red in (a)]. The disappearance of the "red" 
neighbors beside the least fit site by the BS-mechanism 
shifts the overall system to green (more cooperative) with 
time. 

For a quantitative analysis, we measure the mean CP 
(MCP), C(t) = (^Ei c i(*)>i of the populations and 
show the results in Fig. 2. Here, (•) represents the ensem- 
ble average over many different realizations of random 
initial configurations. Note that the MCP also represents 
the overall fitness F(t) = (jj J2i /*(*)) °f the population 



because it is linearly related to MCP: 

F W = ^(E 5 ^-iW + c l+1 (t)] + 2[l- Cl (t)]) 

i 

= 2 + 2(b-l)C(t). (2) 

In Fig. [2j the MCPs for four different system sizes, 
N = 32, 64, 126, and 256, are shown for two different 
values of ui, 0.2, and 1. We use b = 1.5 for all fig- 
ures in this paper, and all data are obtained from nu- 
merical simulations. Because we have assigned a ran- 
dom CP initially, the MCP starts from 0.5 at t = 0. 
For uj = 0.2, the MCP decreases at the beginning and 
then increases to the steady values while it monotoni- 
cally increases from the beginning for uj = 1, as shown 
in Figs. 2(c) and (d). Note that we have two different 
elements in MCP changes. Replacement of the least fit 
member (which is likely to have a high CP) tends to cause 
the MCP to decrease while the replacement of its neigh- 
bors (which probably have low CPs) likely results in an 
increased MCP. The competition between these two el- 
ements governs the early dynamics of the MCP. It can 
decrease initially when uj < uj c « -r = 1/2, where k is 
the number of neighbors. For a sufficiently large sys- 
tem, there would be a site, m, whose CP, c m , is arbitrar- 
ily close to one while those of its neighbors, c m ±i, are 
almost zero. Hence, the expectation of MCP changes, 
AM CP would be ± [(± - l) + kui (± - 0)] and becomes 
negative for uj < r at the beginning. However, as time 
proceeds, the CP values develop spatio-temporal corre- 
lations, and they govern the long-time dynamics. Ini- 
tially, the isolated high-CP cooperators are likely to be 
the least fit member, and they are removed as time pro- 
ceeds. Then, surviving cooperators remain in the groups, 
and, thus, have high fitness. Now, low-CP defectors can 
be the least fit member, especially when they are next to 
a very low-CP member. The replacement of these low- 
CP member by new members with random CPs causes 
the MCP to increase. Therefore, at a later time, the 
MCP easily becomes larger than the initial 0.5 even for 
uj < ui c . Now, a new comer with a random CP will have 
a lower CP than the increased average of MCP. This in 
turn causes the least fit site to be likely located next to 
the new comer's site, resulting in avalanches of replace- 
ment activities. 



IV. ANALYSIS OF THE INITIAL DYNAMICS 

We start from the population with random strategies. 
Hence, there is no correlation between the CPs initially, 
and we may understand the initial dynamics through the 
mean-field calculation. We first define the mean CP of 
the replacement sites (before the replacement), 

~{pmin + 2wC ne i), (3) 



1 + 2ui 



where mean-field dynamics can be easily analyzed. Here, 
Cmin is the average of the CPs for the least fit members, 
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FIG. 3: Evolution of the CPs of the least fit members, Cmin, 
their neighbors, C„ e i, and members that are replaced, C rep , 
are shown together with MCP, C, for (a,c) uj = 0.2, and (b,d) 
uj = 1.0. In (a) and (b), the initial transient behaviors are 
shown while overall behaviors are shown in (c) and (d). The 
system size N = 128 is used for all cases. 



and C ne i is that for the neighbors of the least fit mem- 
bers. On average, CPs of 1 + 2uj sites are updated each 
time. Since the average of the newly assigned random 
cooperation rate is 0.5, C rep satisfies, 



dC, 



rep 



i 



dt 



1 + 2uj 

o.5 - a 



[(0.5 - C mm ) + 2w(0.5 - C nei )] 



rep • 



(4) 



We measure C rep and present them in Fig. [31 together 

tha 



with the CPs of the least fit members, C m i n , that of the 
replaced members, C rep , and the MCP, C, for uj = 0.2 
and uj = 1.0. The C rep curves are, indeed, well described 
by Eq. Q. If we represent the numerical solutions of 
Eq. ([4]) in the figure, they cannot be distinguished from 

th.6 rep 

curves from the simulations because they are 
almost identical. From Fig. [31 we also see that C rep en- 
ters its steady value in a relatively short period of time 
compared to C and rapidly converges to its steady-state 
value of 0.5. For uj — 0.2, the initial C rep is more than 
half and hence decreases to the steady value of 0.5 while 
it increases from the value below 0.5 for u> = 1.0. For a 
sufficiently large system, the initial value of C m i n would 
be 1 while C ne i is 0. Hence, the initial value of C rep would 
be T~rn — 1 which is more than 0.5 for to < 1/2. In this tran- 
sient time of C rep , the dynamics of MCP, C, would be 
mainly determined by the dynamics of C rep . Therefore, 
C initially decreases for uj < 1/2 as does C rep . How- 
ever, after C rep reaches a steady value, the correlation 
of the replacement sites mainly governs the dynamics, 
and C begins to increase. Let m be the least fit mem- 
ber at time t — 1; then, at time t, c m is always updated, 
and Cj =m ±i are updated with probability uj. After re- 
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FIG. 4: Fitness distribution d(f) in the steady states for five 
different system sizes of iV = 16, 32, 64, 128 and 256 with 
(a) uj — 0.2 and (b) uj = 1.0. The system size dependence of 
the effective lower and upper thresholds /£ and (defined 
in the text) are shown in (c). Legends of (a) are also applied 
to (b). 



placement, if the sum of the CPs at these three sites, 

Srep{t) = C m (t) + J2 3 =m±l C j (*) ( at tne tmie ^ S sman > 

at least one of m — 1, morm+1 sites, is likely to have 
small fitness. Therefore, they will be easily replaced in a 
relatively short time. In other words, a new born member 
with small s rep (t) has a short lifetime and contributes less 
to the C than those with large s rep (t). This mechanism 
makes C increase up to (almost) C m i n , and hence, the 
system becomes cooperative overall. Thus, according to 
our model, the emergence of cooperation is intrinsically 
related to the dynamics leading to self-organized critical- 
ity (SOC). 



V. SELF ORGANIZED CRITIC ALITY 

We now show that our model, in fact, drives the popu- 
lation into a SOC state as in the original BS model. We 
measure the distributions of avalanche sizes and distances 
between successive least fit sites in the steady states and 
show that they follow power-law distributions. 

Following Bak and Sneppen [TU], we would like to de- 
fine the size of an avalanche as the number of subsequent 
replacements at the least fit sites below the lower thresh- 
old Jl in its fitness value. The fitness distributions d(f) 
share some characteristics of the BS model [TO] although 
their overall shapes are quite different. A crucial similar- 
ity is that the fitness distribution d(f) in the steady state 
becomes zero for fitness / smaller than a lower threshold 
/l as the system size goes to infinity. 



5 



. 1 1 II 1 

-\ 

r 




(a): 


■ f L =o 


41 + ^ 




r f L=0 


42 




. f L=0 
. . .1 . 


43 * 




1 10 


10 2 ? 10 3 


10 4 








10"' 



d(s) 



10" 



10- 



FIG. 5: (a) Distributions d(s) of avalanche sizes s. Distribu- 
tions with three different values of }l, }l = 2.41, /l = 2.42, 
and /l = 2.43 are measured in systems of N = 256 in their 
steady states. Data with fz, — 2.42 show a most persis- 
tent straight line in the log-log scale fit, indicating the lower 
threshold f L = 2.42 for the N = 256 system with u = 1. The 
black line is the least-squares fit of the data for fz — 2.42 
and is given in a form of d(s) ~ s~ T with r = 0.89 ± 0.05. 
(b) A distribution d(x) of the distances between successive 
minimum fitness sites in the steady states for the system of 
N = 2048 with tj — 1. The black line is the least squares fit 
of the data in the form of d(x) = ax~ a with a = 3.17 ± 0.03. 



The fitness distributions in the steady states for five 
different system sizes are shown in Figs. j4^a) and (b) for 
uj = 0.2 and ui = 1.0. As the system sizes increase, the 
peak positions of the fitness distribution move to the right 
to high values, and the peak widths become narrow. To 
estimate the threshold values ft and fu, we define the 
effective lower [upper] threshold f£(N) [fu(N)] as the / 
value below [above] which the integrated distribution is 
5 percent. We plot them against 1/N in Fig. ^c) for two 
different chain-death rates, ui = 0.2 and uj — 1.0. There 
are no noticeable differences in the thermodynamic values 
for the two uj values. Using linear fitting, we get rough 
estimates of the threshold values, fu — 2.4 ± 0.05 and 
fu = 3.1 ± 0.1, for both u> values. 

For the avalanche size distribution d(s), we need a 
more precise value of //> We measure d(s) with several 
different values of ft around the estimated value. If the 
system is really in a SOC state, we expect the avalanche 
size distribution d(s) to show a power-law distribution, 
for the exact value of ft for the given system. Figure[5f[a) 
shows the distribution of avalanche sizes in a system of 
size N = 256. We plot d(s) against s on a log-log scale 
with three different values of Jl around the value esti- 
mated from Fig. |4|c) to pinpoint the threshold ft . For 
uj = 1.0 shown in Fig. [5ja), the avalanche size distribu- 
tion is well fit by a power-law with ft = 2.42. It remains 
as a line in the log-log plot up to an avalanche size about 
20000, indicating power-law distributions d(s) ~ s~ r . 
The exponent obtained from a least-square fit of the form 
d(s) = Ay~ T is a — 0.89 ± 0.05. This value is consistent 
with the known exponent of the ID BS model 10J . The 
power law indicates that the evolution occurs in a dy- 
namical criticality (TDl [3T]. We measure the avalanche 



distributions for other uj and b and found the critical ex- 
ponent r to be independent of the benefit-to-cost ratio b 
or the chain-death probability uj. 

We also measured the distance distribution between 
successive least fit sites. Denoting the distance between 
successive minimum fitness sites by y, we plot d{y) in 
Fig. [BJb) . The distance distribution is measured in the 
steady states for the system of N = 2048 with w = 1. 
When the distribution d(y) is plotted against y on a log- 
log scale, it also becomes a line, indicating power-law 
distributions d(y) ~ y~ a with the slop a = 3.17 ± 0.03. 
This exponent is also consistent with the known exponent 
of the ID BS model [TU]. It is notable that our model 
belongs to the same universality class as the BS model 
in spite of the complexity in computing the fitness of 
members and the non-trivial dynamics of the population- 
fitness changes. 



VI. CONCLUDING REMARKS 

We have considered the BS mechanism as a repro- 
duction process with fitness given by a PD game pay- 
off on a network structure. Our observation may have 
more natural implication in economical systems because 
the BS process with chain bankruptcy is a more feasi- 
ble scenario. It might be worthwhile analyzing weekly or 
monthly bankruptcy data and see if they follow a power- 
law distribution as our study suggests. 

We have simulated our model with other values of the 
benefit-to-cost ratio b and see that cooperation emerges 
in a wide range of chain-death rates w, as long as b is 
larger than 1. In contrast to a common belief, cooper- 
ation can emerge even with parameters that a popula- 
tion with random strategies decreases cooperation. This 
is possible because the BS mechanism builds dynami- 
cal correlations that suppress the long-term survival of 
non-cooperators even in the region where mean-field cal- 
culation predicts a decrease in cooperators. The same 
dynamical correlation leads to SOC in replacement activ- 
ities with the same exponents as the original BS model. 
The strategy space presented here is rather small. Mixed 
but only history independent strategies are considered on 
a very simple population structure, a ID lattice. How- 
ever, we speculate that our main results, the emergence of 
cooperation and SOC, are robust under variations in the 
population structure or the strategy space extension. In 
fact, the preliminary results with the extended strategy 
space show that the emergence of cooperation appears 
more easily and rapidly when the reactive strategies are 
included. 



VII. ACKNOWLEDGEMENTS 

This work was supported by the National Research 
Foundation of Korea Grant funded by the Korean Gov- 
ernment (MEST) (NRF-2010-0022474). H.-C. J. would 



6 



like to thank KIAS for the support during the visit. 



[1] M. A. Nowak, Evolutionary Dynamics (The Belknap 

Press of Harvard University Press, Cambridge, 2006) . 
[2] M. A. Nowak, Science 314, 1560 (2006). 
[3] M. Milinski, Nature 325, 433 (1987). 
[4] S. VanSegbroeck, F. C. Santos, T. Lenaerts, and J. M. 

Pacheco, Phys. Rev. Lett. 102, 058105 (2009). 
[5] J. Gomez-Gardenes, M. Campillo, L. M. Floria, and Y. 

Moreno, Phys. Rev. Lett. 98, 108103 (2007). 
[6] J. M. Pacheco, A. Traulsen, and M. A. Nowak, Phys. 

Rev. Lett. 97, 258103 (2006). 
[7] F. C. Santos and J. M. Pacheco, Phys. Rev. Lett. 95, 

098104 (2005). 

[8] S. J. Gould and N. Eldredge, Paleobiology 3, 115 (1977). 

[9] D. M. Raup, Science 231, 1528 (1986). 
[10] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993). 
[11] W. D. Hamilton, J. Theor. Biol. 7, 1 (1964). 
[12] M. A. Nowak and K. Sigmund, Nature 393, 573 (1998). 
[13] S. A. West, I. Pen, and A. S. Griffin, Science 296, 72 



(2002). 

[14] A. Grafen, J. Evol. Bio. 20, 2278 (2007). 

[15] J. A. Scheinkman and M. Woodford, The American Eco- 
nomic Review 84, 417 (1994). 

[16] R. V. Sole and S. C. Manrubia, J. Theor. Biol. 173, 31 
(1995). 

[17] T. Killingback and M. Doebeli, J. Theor. Biol. 191, 335 
(1998). 

[18] A. Arenas, A. Diat-Guilera, C. J. Perez, and F. Vega- 
Redondo, J. Econ. Dyna. & Cont. 26, 2115 (2002). 

[19] H. Ebel and S. Bornholdt, Phys. Rev. E 66, 056118 
(2002). 

[20] J. A. Fletcher and M. Doebeli, Proc. R. Soc. B 276, 13 
(2009). 

[21] P. Bak, How Nature Works (Springer- Verlag, New York, 
1996). 



